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I.  INTRODUCTION 


Electromagnetic  wave  propagation  over  hilly  terrain  is  important  not  only  in  point-to- 
point  communication  over  land,  but  also  in  ground-to-air  communication.  Of  late,  it  has 
assumed  importance  in  outdoor  propagation  in  the  context  of  personal  communications 
network  design.  An  exact  analytical  solution  of  the  problem  for  general  terrain  is  not  pos¬ 
sible  and  one  often  resorts  to  an  approximate  or  a  numerical  approach.  In  a  previous  paper 
[1],  we  developed  a  numerical  model  for  propagation  predictions  over  inhomogeneous,  ir¬ 
regular  terrain  using  the  magnetic  field  integral  equation.  Although  the  method  could 
include  all  effects  of  wave  propagation  such  as  reflection,  diffraction,  surface  wave  excita¬ 
tion,  and  backscattering,  a  principal  limitation  of  the  method  was  the  requirement  of  large 
computer  resources  (CPU  time  and  memory),  particularly  for  electrically  large  terrain  ir¬ 
regularities.  For  instance,  if  the  integral  equation  is  solved  numerically  by  the  method  of 
moments  [2],  the  matrix  fill  time  would  be  of  order  0(N'^)  and  the  inversion  time  of  order 
0{N^),  where  N  is  the  total  number  of  unknowns.  As  the  matrix  generated  is  dense,  the 
memory  requirements  would  be  of  order  0(N‘^).  The  method  is  attractive  for  small  terrain 
irregularities,  but  is  computationally  prohibitive  for  large  terrain  irregularities  [1]. 

In  this  paper  we  present  a  computationally  efficient  model  of  wave  propagation  based 
on  finite  differences.  It  may  be  noted  that  this  method  has  no  semblance  to  the  one 
proposed  in  [3],  where  one  proceeds  with  the  parabolic  equation  approximation  of  the 
Helmholtz  equation.  Use  of  parabolic  equation  approximation  precludes  backscattering, 
which  is  sometimes  important.  In  contrast,  we  apply  finite  differences  directly  to  the 
Helmholtz  equation  without  introducing  any  dubious  approximations.  Even  though  one 
has  to  handle  a  larger  matrix  when  dealing  with  finite  differences  in  contrast  to  boundary 
methods  such  as  integral  equation  methods,  the  resulting  matrix  is  sparse,  and  often  faster 
to  invert  than  the  dense  matrix  generated  via  the  latter.  The  matrix  fill  time  is  still  of 
order  O(N^),  but  the  inversion  time  is  reduced  to  a  lower  order  of  0{N).  Thus,  the  real 
advantage  of  a  finite  difference  scheme  is  felt  for  large  problems  where  the  run  time  is 
dominated  by  the  inversion  time.  Substantial  savings  in  memory  are  also  be  accomplished 
as  the  matrix  is  sparse.  The  memory  required  is  only  of  order  0{N).  The  method  can 
potentially  solve  larger  problems  than  possible  with  the  method  of  moments.  The  principal 
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difficulty  associated  with  the  use  of  finite  differences  is  the  treatment  of  mesh  truncation 
when  applied  to  open  type  of  problems  such  as  encountered  in  propagation,  antennas,  and 
scattering.  We  will  use  a  method  known  as  the  Measured  Equation  of  Invariance  (MEI) 
originally  proposed  by  Mei,  et  ah,  [4]  to  treat  the  mesh  truncation. 

The  term  scatterer  will  be  used  generically  to  denote  an  obstacle  in  an  open  envi¬ 
ronment.  Local  radiation  conditions  such  as  those  proposed  in  [5]  and  [6]  are  good  when 
the  truncating  boundary  is  far  removed  so  that  it  is  in  the  far-zone  field.  However,  this 
results  in  a  larger  matrix  size  with  obvious  implications  to  computational  time  and  mem¬ 
ory.  Global  type  radiation  conditions  [7]  can  be  used  on  a  tighter  terminating  boundary 
to  result  in  a  smaller  computational  domain.  However,  this  will  destroy  the  sparsity  of  the 
matrix  generated  by  finite  differences  and  defeats  the  whole  purpose  of  employing  it  in  the 
first  place.  What  is  needed  is  a  local  boundary  condition  of  the  type  in  [5],  but  applicable 
very  close  to  the  scatterer.  Although  it  is  far  more  complex  to  find  near  field  radiation 
conditions  than  it  is  to  find  far  field  ones,  it  is  partially  offset  by  the  fact  that  boundary 
conditions  on  a  continuous  spatial  domain  are  not  needed  when  finite  methods  are  used. 
The  MEI  method  enables  one  to  generate  the  near-field  conditions  over  a  discrete  domain. 

In  this  paper  we  deal  with  only  two-dimensional  sources  and  fields,  and  one  dimen¬ 
sional  terrain  characteristics.  Accordingly,  the  terrain  properties,  the  sources,  and  the 
corresponding  fields  are  all  invariant  with  respect  to  the  longitudinal  variable  z.  It  is  as¬ 
sumed  that  the  terrain  is  characterized  by  its  local  impedance  and  height  over  a  reference 
plane,  both  of  which  may  vary  from  point  to  point.  In  section  Ila,  we  present  a  finite  dif¬ 
ference  discretization  of  the  two  dimensional  Helmholtz  equation  and  present  an  overview 
of  the  present  method.  To  realize  the  mesh  termination  conditions  via  the  MEI  method 
described  in  Ild,  an  accurate  representation  of  the  near-zone  scattered  field  is  necessary. 
In  section  Ilb  we  give  an  integral  representation  of  the  scattered  field  of  a  TE^  line  somce 
over  an  irregular,  impedance  surface.  The  corresponding  expressions  for  the  incident  field 
and  the  Green’s  functions  are  presented  in  section  lie.  In  section  He,  we  describe  the  Ri- 
catti  block-by-block  elimination  technique  [8]  of  sparse  matrix  inversion.  Finally  in  section 
III  we  present  numerical  results  for  both  the  sky  wave  and  the  ground  wave  and  provide 
comparisons  for  test  geometries. 
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II.  FORMULATION 


In  the  present  paper,  we  consider  a  two  dimensional  situation  as  shown  in  Fig.  1.  The 
transmitting  antenna  is  a  transversely  polarized  electric  line  source  located  at  {xotVo)- 
Such  a  source  will  have  its  electric  field  confined  to  the  transverse  xy  plane,  and  is  the 
two-dimensional  counterpart  of  the  superposition  of  a  vertical  electric  dipole  (VED)  and 
a  horizontal  electric  dipole  (HED)  in  three  dimensions.  The  field  due  to  the  source  can 
be  classified  as  TE^,  and  all  components  could  be  expressed  in  terms  of  the  ^-component, 
of  the  magnetic  field.  It  is  assumed  that  all  distance  variables  are  normalized  with 
respect  to  the  free-space  wavenumber  kg  =  where  u:  is  the  radian  frequency  of 

the  wave,  Cq  is  the  permittivity,  and  fig  the  permeability  of  free-space.  Accordingly,  we  set 
X  :=  kgX,  y  kgy,  etc.  An  time  dependence  is  assumed  and  suppressed.  For  TE^ 
polarization,  the  impedance  boundary  condition  [9]  of  the  form  h  x  E  =  rjoAh  X  (n  x  H), 
relating  the  electric  field  vector  E  to  the  magnetic  field  vector  H  leads  to 


dHz 

dn 


(1) 


where  the  unit  normal  n  points  out  of  the  impedance  surface,  r)o  =  yiXg/tg  is  the  intrinsic 
impedance  of  free  space,  and  A  is  the  normalized  impedance  of  the  surface. 


Ila.  Overview  of  the  Method: 

We  let  if  to  denote  the  ^-component  of  the  scattered  magnetic  field  due  to  a  TE^  line 
source  over  an  inhomogeneous,  irregular  terrain.  We  assume  that  if  represents  scattering 
only  from  the  irregularites  and  inhomogeneities  in  a  reference  impedance  plane.  Thus, 
the  scattered  field  is  identically  zero  when  the  terrain  is  fiat  having  an  impedance  equal 
to  the  reference  impedance.  The  scatterer  then  consists  of  those  portions  of  the  terrain 
where  (i)  the  impedance  is  different  from  the  reference  value,  and  (ii)  the  elevation  is 
different  from  zero.  The  computational  domain  consists  of  a  region  in  space  bounded  by 
the  terrain  at  the  bottom  and  a  terminating  (Le.,  truncating)  boundary  at  the  top  as 
shown  in  Fig.  2.  It  is  assumed  that  the  boundary  of  the  scatterer  is  subdivided  into  N  —  1 
segments,  thereby,  generating  N  points  on  it.  The  terminating  boundary  is  similarly 
partitioned  into  —  1  segments.  We  generate  a  structured  mesh  in  the  computational 
domain  by  adding  M  interior  layers  between  the  object  boundary  and  the  terminating 
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boundary,  each  having  N  points.  A  total  of  M  +  2  layers,  each  having  N  points  are 
thus  generated.  Layer  numbering  is  done  in  ascending  order  starting  from  the  object 
boundary  (layer  number  0)  and  progressing  towards  the  terminating  boundary.  Node 
numbering  is  done  from  left  (number  1)  to  right  (number  N).  We  use  the  notation 
m  =  0,l...,M  +  l,  n  =  l,...,iVto  denote  the  cartesian  coordinates  of  the  nth  point  on 
layer  number  m.  Points  on  the  layer  immediately  following  the  scatterer  are  assumed  to 
lie  on  the  local  normal  emanating  from  a  point  on  the  scatterer.  This  is  to  permit  easy 
implementation  of  the  impedance  boundary  condition  of  (1).  Within  the  computational 
domain,  the  scattered  field  i/>  satisfies  the  scalar  Helmholtz  equation 


a2  \ 


Since  the  terrain  irregularities  do  not  necessaily  conform  to  any  one  standard  coordinate 
system,  the  mesh  is  non-orthogonal,  and  traditional  finite  difference  equations  are  not 
applicable.  We  will  develop  the  required  finite  difference  equations  similar  to  the  method 
outlined  in  [10].  Prompted  by  the  presence  of  second  order  derivatives  in  the  Helmholtz 
equation,  we  choose  a  five-point  star  mesh  centered  at  0(Xo,yo)  and  sorrounded  by  four 
neighboring  nodes  with  local  coordinates  (Xfc,^),  ^  =  1, ...  ,4,  as  shown  in  Fig.  3b.  The 
global  indices  of  the  nodes  0,  1,  2,  3,  and  4  are  (m,  n),  (m  -  1,  n),  (m  -f  1,  n),  (m,  n  -  1),  and 
(m,n  -|-  1)  respectively.  Using  the  notation  in  the  local  coordinates  that  =  rp{Xk,Yk), 
we  assume  a  finite  difference  equation  over  the  mesh  in  the  form 


4 

V’o  +  =  0, 

Jt=l 


(3) 


where,  Ck,  are  unknown  complex  constants.  The  above  equation  can  be  rewritten  as 


Y^Ck 


k=l 


th 

V’o 


=  -1. 


(4) 


The  coefficients  are  determined  by  choosing  four  linearly  independent  plane  waves  (i.e., 
trial  solutions  of  (2))  traveling  along  the  lines  joining  the  central  node  to  the  four  neigh¬ 


boring  nodes;  i.e.,  we  choose 


^  _  g-jd*  cos(ait-Q!,)  4  (5) 

V’o 
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where  (dk,ak)  are  the  polar  coordinates  of  the  neighboring  nodes  with  respect  to  the 
central  node  0.  The  linear  system  of  equations  resulting  from  (4)  can  then  be  solved 
for  the  unknowns  c^.  It  is  interesting  to  note  that  the  values  for  the  coefficients  reduce 
to  the  standard  values  when  the  mesh  conforms  to  a  standard  coordinate  system.  The 
computational  domain  is  open  ended  at  the  sides  and  one  cannot  choose  a  five-point  star 
mesh  as  above.  At  the  two  ends  of  the  domain,  a  modified  mesh  as  shown  in  Figs.  3a  and 
3c  is  used.  For  the  mesh  shown  in  Fig.  3a,  it  is  possible  that  the  plane  waves  traveling 
along  directions  3  and  4  become  linearly  dependent,  or  nearly  so.  To  avoid  this  degeneracy, 
we  reverse  the  direction  of  the  plane  wave  traveling  along  0  —  3  and  solve  the  linear  system. 
Likewise,  to  avoid  degeneracy  with  the  mesh  of  Fig.  3c,  we  reverse  the  direction  of  the 
plane  wave  traveling  along  0  —  4  and  solve  for  the  coefficients. 

At  the  lower  (object)  boundary  of  the  computational  domain,  the  impedance  boundary 
condition  (1)  is  applicable  to  the  total  field  =  X  +  where  x  =  is  the  incident 
field.  For  small  distances,  the  solution  of  (1)  along  the  normal  results  in 

«  =  -XS  +  =  -Xo”  +  («  +  Xf  )/>.,  (6) 

where  hn  =  -f  (yf  -  yj the  normal  distance  between  scatterer  and  the 

first  interior  layer  at  the  nth  node.  A”  is  the  value  of  at  the  nth  node,  Vm), 

Xm  = 

Additional  conditions  are  needed  at  the  artificial,  terminating  boundary.  To  simulate 
free-space,  boundary  conditions  proposed  by  Mei  et  al.,[4]  are  assumed  to  hold  on  the 
boundary.  Accordingly,  on  the  five-point  star  mesh  of  Fig.  3b  formed  by  layer  numbers 
M  —  1,  M,  and  M  +  1,  the  scattered  field  satisfies  (in  the  local  coordinates) 

4 

V'o  +  ^  ajfeV’fc  =  0,  (7) 

k=l 

where  the  complex  coefficients  axe  different  from  Cfc.  They  are  determined  in  a  fashion 
similar  to  the  latter  except  that  the  trial  solutions  are  no  longer  plane  waves  but  are 
dependent  on  the  geometry  of  the  scatterer  and  the  location  of  the  mesh.  Such  trial 
solutions  can  be  obtained  from  the  near-zone  behavior  of  the  field  and  further  discussed 
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in  section  Ild.  It  is  important  to  note  that  a  discrete  equation  of  the  form  (7)  over  a  five- 
point  computational  molecule  of  Fig.  3  can  only  capture  spatial  derivatives  upto  second 
order  (without  the  cross  terms).  Equations  (3),  (6)  and  (7)  are  combined  to  result  in  a 
matrix  equation  of  the  form 

-f  -1- (8) 

+  =F",  2<n<N-l  (9) 

+Bn^^  +  =  (10) 


In  the  above  matrices,  c^’”  is  the  kih.  finite  difference  coefficient  associated  with  the  node 
at  {x^,y^),  and  is  the  A:th  MEI  coefficient  associated  with  the  node  at  ).  The 

system  of  equations  defined  by  (8—10)  can  be  solved  efiiciently  by  the  Ricatti  block-by- 
block  elimination  technique  [8]  described  in  lie.  After  the  scattered  field  is  solved  in  the 
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interior  domain,  the  total  field  on  the  scatterer  can  be  recovered  from  (6).  The  field  at  any 
point  can  then  be  determined  from  the  total  field  on  the  scatterer  by  using  the  integral 
representation  given  in  the  next  section.  In  the  next  few  sections,  we  give  details  of  the 
various  steps  presented  in  this  section. 

Ilb.  Integral  Representation: 

Fig.  1.  shows  an  electric  line  source  Jq  located  at  0(xo,yo)-  If  is  required  to  find 
the  total  field  {E,H),  at  an  arbitrary  point  F{xp,yp)  over  the  irregular,  inhomogeneous 
terrain.  We  will  employ  the  reciprocity  theorem  [2]  to  express  the  total  field  in  terms  of 
the  scattered  field  and  an  incident  field,  (Eo.Ho),  defined  to  exist  over  a  flat,  constant 
impedance  plane,  Cq,  of  normalized  impedance  Ao-  The  general  relationship  between  a 
surface  impedance  A  and  the  local  ground  constants  (coCr,  A'o,  o’)  is 

A  =  l/\/(er  (15) 

where  cTr  is  the  relative  conductivity  representing  <r/a;eo)  o’  being  the  actual  conductivity 
of  earth.  The  above  relationship  is  used  to  evaluate  not  only  Aq,  but  also  A^  that  will 
be  encountered  below.  The  field  {E,H)  is  different  from  {Eo,Ho)  due  to  (a)  the  terrain 
inhomogeneities,  or  (b)  terrain  irregularities,  or  (c)  both.  Let  {Ei,Hi)  be  the  fields  due  to 
sources  (Ji,Mi)  located  at  F{Xp,yp)  over  the  reference  impedance  plane.  The  fields  ‘0’  and 
T’  satisfy  the  impedance  boundary  condition  y  x  Eo^  =  rjoAgy  x  {y  x  Ho)  on  Cq,  where 
y  is  the  unit  normal  on  Cq.  The  total  fields  satisfy  the  impedance  boundary  condition 
y  X  E  =  rio^s'iT'  X.  (h  X  H)  on  the  terrain  surface  C,  where  the  impedance  A^  is  possibly 
different  from  Ao,  and  may  vary  from  point  to  point.  The  surface  C  is  assumed  to  consist 
of  flat  portions  Cj  where  the  impedance  A^  deviates  from  Ao  and/or  Ct  where  the  terrain 
deviates  from  being  flat.  The  remaining  portion  Cr  of  the  surface  is  concident  with  the 
corresponding  portions  of  Cq.  The  surface  Ct  is  the  projection  of  Ct  onto  Cq.  Let  Coo 
be  a  semicircular  cylinder  bounded  at  infinity,  and  A  be  the  interior  region  bounded  by 
the  surfaces  C  and  Coo.  The  unit  normal  n  points  into  the  region  A.  It  is  assmned  that 
the  ground  deviates  from  the  reference  plane  only  in  the  positive  y  direction.  Thus,  we 
consider  only  hilly  obstacles  and  do  not  address  the  problem  of  trenches  in  a  flat  plane. 
The  ground  constants  are,  however,  allowed  to  have  variations  along  the  terrain. 
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The  fields  {E  -  Eo,H  -  Ho)  have  no  sources  in  the  region  A,  and  satisfy  source-free 
Maxwell’s  equations,  while  the  fields  i^E\^H\)  satisfy  Maxwells  equations  with  sources 
From  a  two  dimensional  version  of  the  reciprocity  theorem  [2],  we  have 


J  ^{E  -  Eo)  X  Hi  -  El  X  {H  -  Ho)  •nd£  =  Jj [(E  -  ^o)  ■  Ji  -  (H  -  Ho)  ■ 


Ml  da 


C’r+C'^+Ct+Coo 


The  contour  integral  over  Coo  vanishes  due  to  the  Sommerfeld  radiation  conditions  [11]. 
Since  (E,H)  and  {Eo,Ho)  satify  the  same  boundary  conditions  on  Cr,  the  integral  over  it 
vanishes.  Applying  reciprocity  theorem  to  the  source-free  region  bounded  by  Ct  —  Cb,  we 


obtain 


Eo  X  Hi  -  El  X  Ho  -hdi-O 


Ct-Ct, 


Hence 


Eo  X  Hi  -  El  X  Hq  •  n 


ndi  =  J  j^.Fo 


X  Hi  -  El  X  Ho  •  n  =  0. 


The  last  equality  follows  from  the  fact  that  the  ‘0’  and  T’  fields  satisfy  the  same  impedance 
boundary  conditions  on  Cj.  Making  use  of  these  in  (16)  we  get 

ff  \{E-Eo)- Ji-{H- Ho) -Ml  da=  f  {E  x  Hi  -  Ei  x  H)  ■  h  dt  (19) 

A  Cj+C, 

Next,  we  pick  a  f-directed  magnetic  line  source  of  unit  voltage  for  Mi  and  choose  Ji  =  0 


to  arrive  at 


=  [(li  X  B)  ■#!'>  + (ft  X  5)  ■£''>]  <«,  (20) 

where  {e[''\  H^'"'^)  is  the  field  due  to  a  z-directed  magnetic  line  source  of  unit  voltage 
located  at  the  field  point  P.  We  relate  (h  x  H)  and  (n  x  E)  on  Cf  -i-  Ct  by  means  of  the 
impedance  boundary  condition  to  finally  obtain 

z- H{p)  =  z- Ho{p)-J  E["\Q)-r}oAs(nx  h["\q)^Y  (nx  H(q)'^  di^.  (21) 
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The  above  equation  is  rewritten  as 


H,{p)  =  Hoz{p)  -  j  H,{Q)g{Q,p)  di^ ,  (22) 

C=Cj+Ct 

where 

g{Q,p)  =  E[f{Q)  -  r]oAsH[l\Q)  (23) 

is  defined  as  the  Green’s  function  at  Q  due  to  a  z-directed  magnetic  line  source  located  at 
P  over  a  plane  of  impedance  Aq.  The  unit  tangent  1.  to  the  contour  is  defined  such  that 
i  X  h  =  z.  Equation  (22)  may  be  used  to  set  up  an  integral  equation  for  the  unknown 
Hz  and  solved  numerically.  However  this  will  be  computationally  intensive  as  already 
discussed  and  will  not  be  considered  in  the  present  report.  In  the  next  section,  we  present 
expressions  for  Hqz  and  g{Q,P)  which  allow  rapid  numerical  computation. 

lie.  Incident  Magnetic  Field  and  Green’s  Function: 

Fig.  4.  shows  a  transversely  polarized  electric  line  source  located  at  (xo,yo)  over 
a  plane  of  constant  impedance  Aq.  Due  to  the  impedance  boundary  condition  we  have 
Eox  =  tjoAoHqz-  The  line  source  carries  a  total  current  Ig  and  has  its  current  moment 
directed  along  the  unit  vector  i  which  makes  an  angle  9o  with  the  x-axis.  The  current 
moment  corresponding  to  the  mirror  image  of  the  source  about  a  perfectly  conducting 
plane  at  y  =  0  points  in  a  direction  i'  which  makes  an  angle  ir  —  6o  with  the  x-axis.  The 
current  density  of  the  line  source  is  expressed  as 

J  =  iIoklS{x  -  Xo)6{y  -  yo)  kll,  (24) 

where  ^(•)  represents  a  unit  delta  function.  The  magnetic  vector  potential  A  =  xAx  +  yAy 
satisfies 

V^A  +  A  =  -f^oI.  (25) 


It  is  easily  seen  from  Maxwell’s  equations  that 

'dAy  dAx 
dx  dy 
dHox 


UJ 


Hez  =  - 
yo 

Eqx  =  -jvo 


dy 


(26) 

(27) 
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Equation  (25)  can  be  seperated  into  its  respective  cartesian  components.  A  vertical  current 
excites  a  vertical  component,  Ay,  of  the  vector  potential,  whereas,  a  horizontal  component 
excites  only  a  horizontal  component,  A^.  It  is  to  be  noted  here  that  a  horizontally  polarized 
line  source  will  only  need  for  the  satisfaction  of  the  boundary  conditions  at  y  =  0.  This 
is  in  contrast  to  the  three-dimensional  case  where  a  HED  requires  both  and  A^  for  the 
satisfaction  of  boundary  conditions.  We  take  a  Fourier  transform  on  both  sides  of  (25) 
with  respect  to  the  x-axis,  and  denote  the  transformed  quantities  with  a  tilde  and  the 
normalized  transform  variable  by  k^.  When  this  is  done,  an  inhomogeneous  harmonic 
equation  with  respect  to  y  is  obtained  in  the  transformed  domain.  The  components  of  the 
vector  potential  can  be  determined  from  it  as 


Ax(^i ,  y)  )  yo) 
Ay  (A^i ,  y,  yo) 


HoIo  cos  000^^^^“ 

HqIq  sin 

4:Trjl3 


je-l^ly-s/ol  +  ^  (29) 


where  Rh{/3)  and  Rvi/3)  are  the  reflection  coefficients  for  the  horizontal  and  vertical  polar¬ 
izations  respectively,  and  =  1.  Imposition  of  the  impedance  boundary  condition 

at  y  =  0  yields 

=  =  (30) 

The  components  of  vector  potential  in  the  space  domain  are  obtained  by  taking  the  inverse 
Fourier  tranformation  of  (28)  and  (29).  Letting  the  free-space  Green’s  function  11  as 


OO 

Il{x,y,Xo,yo)  =  J 


,-j[kx(x-xo)+^\y-yo\ 


dkx  =  Ho^^[\/{x  -  XoY  +  (y  -  yo)^],  (31) 


we  arrive  at 


A  =  |zn(2:,y;Xo,yo)  +  i'  [n(x,y;Xo,  -yo)  -  2Ao(5n  (x,  y;  Xo,  yo)]| ,  (32) 


where 


{x,y;xo,yo)=  /  - - 


-3[0(.y+yo)  +  kx(.i:-Xo)] 


/5(/?  +  Ao) 


The  integral  defined  in  (33)  is  of  the  so  called  Sommerfeld  type  and  difficult  to  evaluate 
directly.  To  efficiently  compute  the  incident  magnetic  field  and  the  Green’s  function,  we 
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will  incorporate  the  ideas  contained  in  [12]  and  [13].  For  small  distances,  we  use  a  modified 
method  of  [13].  For  large  horizontal  distances,  where  the  method  of  [13]  fails,  we  have  used 
a  modified  version  of  the  method  in  [12]. 

For  small  distances,  we  employ  the  clever  device  originally  used  by  Chow,  et  al,  [13] 
to  calculate  the  Green’s  function  for  microstrip  structures.  The  key  to  the  procedure  is 
to  convert  (31)  to  an  integral  over  the  complex  /9-plane  and  make  approximations  to  (33). 
Firstly,  we  note  that  the  imaginary  part  of  /3  must  be  chosen  appropriately  (i.e.,  9=(^)  <  0) 
such  that  integral  in  (31)  converges  for  all  jy  -  y^j  >  0.  By  a  change  of  variable,  (31)  may 
be  converted  to  an  integral  in  terms  of  /9  as 

,  ,  _  o  r  f,-}^\y-yo\ 

^?  [\/(^  -  +  {y-  yo)^]  =  -  J  - j, - cos[k^(x  -  Xo)]  d/9,  (34) 

Ti 

where  the  contour  of  integration,  Fi,  is  shown  in  Fig.  5.  Similarly, 

2  f 

6U  (x,  y;  Xo,yo)  =  -  J  ^  ^  ^  cos[k^(x  -  Xo)]  djS.  (35) 

Ti 

The  contour  of  integration  for  <911  may  be  deformed  to  a  modified  contour,  r2,  since  no 
singularities  are  enclosed  between  it  and  Fi  (see  Fig.  5).  The  contour  F2  consists  of  two 
parts:  a  straight  line,  Fg,  drawn  from  (1,0)  to  (0,  —jTo),  for  some  chosen  positive  constant 
To,  and  a  portion  coincident  with  that  of  Fi.  Most  of  the  contribution  to  the  integral 
comes  from  the  former  as  the  integrand  decays  exponentially  on  the  latter  for  sufficiently 
large  y  -I-  yo-  The  trick  is  to  approximate  the  spectral  function  l/(/9  -|-  Aq)  over  Fg  as  a 
finite  sum  of  complex  exponentials  of  the  form 


1 

(^  +  Ao) 


Nt 

jfc=i 


(36) 


The  complex  constants  Ak  and  fjt  can  be  determined  by  the  Prony’s  method  [14].  Substi¬ 
tuting  this  into  (35)  we  have 


dn(a:,y;Xo,yo) 


Nt  2  . 

A:=l  T 


-iP{y+yo-jtk) 


cos[fci(x  —  Xo)]  d/9 
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+ 


—joo 

2  f  cos[A:j;(a:  —  Xq)] 


/ 


TT 

/3=-jTo 
Nt 


Nt 


dp 


=  '£,AtHi^'>(Rt)+~»  [z2£i(z2)  -  e-‘’] 


A:=l 
.  Nt 


ttTo 


-^yAk  [Ei{z2-jhTo)  +  E^(z*  -jhTo)] 

‘TT  f  ^ 


(37) 


k=l 


The  first  term  of  the  last  equality  follows  by  comparison  with  (34).  The  quantity  Rk  is  the 
distance  between  the  observation  point  (x,  y)  and  the  complex  image  at  (xq,  —yo'^-j'tk)  given 
by  R},  =  ^J{x  -  XoY  +  (y  +  2/o  -  j^kY,  with  the  square  root  chosen  such  that  ^{Rk)  >  0, 
and  Z2  —  [{y  +  yo)  +  -  ^o)]To.  The  remaining  integral  in  (37)  over  the  negative 

imaginary  axis  has  been  performed  in  closed  form  in  terms  of  exponential  integrals,  Ei{-) 
[21],  by  using  the  approximations  ^  +  Ag  ~  and  kx  ~  —j^-  The  quantity  611  can 
be  calculated  much  faster  using  (37)  rather  than  (33)  or  (35).  However,  when  y  +  yo  is 
small  and  |x  —  Xq]  is  large,  the  contribution  arising  from  the  tail  of  r2  becomes  important, 
and  the  exponential  approximation  to  the  spectral  function  in  (36)  requires  many  terms. 
Using  only  a  few  terms  leads  to  large  errors  in  the  computation  of  6n.  The  accuracy  can 
be  greatly  improved  at  the  cost  of  computational  time  by  expressing  the  spectral  function 
l/(^  +  Ao)  over  as  a  Laplace  integral  of  the  form 


1 

(/?  +  Ao) 


OO 

J  e-^o^e-^Ut. 

t-O 


(38) 


Such  a  procedure  was  suggested  in  [12].  The  Laplace  transform  of  the  spectral  function  in 
(38)  converges  for  -7r/2  -  arg(Ao)  <  arg(t)  <  7r/2  -  arg(Ao).  Substituting  this  into  (35) 
and  changing  the  order  of  integration,  we  arrive  at 

OO 

6n(x,y;Xo,yo)  =  y  e~^°^Hf\D)dt,  (39) 

<=o 


with  the  complex  distance  D  =  -^(x  —  Xg)^  +  (y  +  yo  —  defined  such  that  3i(I?)  >  0 
for  representing  an  outgoing  wave  and  ^{D)  <  0  for  the  convergence  of  the  integral.  Even 
though  it  is  much  faster  to  compute  6n  using  (39)  than  using  (35),  the  integration  can  be 


12 


slow  for  certain  combinations  of  Ao,  {x  —  Xo),  and  (y  +  yo)-  This  is  because  the  integrand 
decays  rather  slowly  on  the  positive  real  axis  of  the  t-plane.  To  further  aid  fast  evaluation 
of  the  integral  (39),  we  deform  the  original  contour  in  the  complex  t— plane  to  one  over 
which  the  asymptotic  form  of  the  Hankel  function  decays  most  rapidly  (steepest  descent 
path).  This  is  obtained  by  setting  ^D(t)]  =  =  0)]  =  r2  on  the  deformed  contour, 

where  r2  =  y/ix  -  +  (y  +  VoY  is  the  distance  of  the  observation  point  from  the  mirror 

image,  li  t  =  u  +  jv,  the  steepest  descent  path  through  t  =  0  can  be  obtained  as 

u — >[0,oo),  V  := -(y  +  yo)  +  r2^J~—^^^^  (40) 

on  which  the  distance  function  is 


Evaluation  of  the  integral  over  the  steepest  descent  path  permits  rapid  computation 
of  ^n.  Although  this  is  true  uniformly  for  cill  x  —  xo  and  y  +  yo,  it  is  less  efficient  than  (37) 
for  small  values  of  r2.  This  is  due  to  the  presence  of  the  Hankel  function  with  complex 
arguments  in  the  integrand  of  (39),  which,  one  has  to  evaluate  repeatedly  to  perform  the 
integration.  Hankel  functions  need  be  computed  only  for  a  fixed  number  of  arguments  in 
(37)  in  contrast  to  the  several  tens  of  times  needed  when  (39)  is  used.  It  is  very  important 
indeed  to  reduce  the  time  required  to  compute  the  incident  field  and  the  Green’s  function, 
particularly  when  one  has  to  evaluate  them  several  thousands  of  times  as  in  the  present 
method.  This  point  will  be  appreciated  shortly  when  we  discuss  the  MEI  method. 

Combining  (26)  and  (32)  the  incident  magnetic  field  can  be  obtained  as 

Ho.  =  =^|sin(«,  -  «.)«r(n)  +  sinCSo  +  h)ni^\r2) 

g 

+  2Ao(sin  ^0  ^  +  cos 

where  =  ^/{x  —  Xq)^  "H  (y  —  yo)^i  &nd  ^1(^2)  Is  the  positive  angle  made  with  the  x-axis 
by  the  straight  line  joining  the  observation  point  and  the  source  (mirror  image)  point. 

Using  a  similar  analysis,  the  Green’s  function  (defined  through  (23))  at  the  point 
Q(xg,yg)  due  to  a  unit  voltage,  z-directed,  magnetic  line  source  located  at  P(xp,yp)  can 


-)OT(x,y;Xo,yo)  >, 
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be  obtained  as 

Q{q,p)  =  ^|sin(^  -  9{)Hf\rz)  +  sin(^  -  -jAg  J^q^Vs)  + 

4;  L 

d  d  "  ') 

+2Ao  jAs  +  sin^^ - cos0—  8Il{xg,yg]Xp,yp)\,  (43) 

OXq  ^Vq  J  J 

where  =  y/{xq  -  XpY  +  (j/,  -  ypY,  r4  —  ^{xg  -  XpY  +  (yg  +  ypY,  ^3(^4)  is  the  positive 
angle  made  with  the  x-axis  by  the  straight  line  joining  the  observation  point  Q  and  the 
source  point  P  (mirror  image),  and  9  is  the  angle  made  with  the  positive  x-axis  by  the 
unit  tangent  at  Q. 

An  integral  equation  may  be  derived  for  the  unknown,  Hz,  on  the  scatterer  by  substi¬ 
tuting  the  above  Green’s  function  into  (22)  and  taking  the  limit  as  the  point  P  approaches 
the  scattering  boundary  from  outside.  When  this  is  done,  one  gets 


Hz{p)  =  2Hoz{p)-2  [  HziQ)g(Q,p)diQ,  P,QeCj+c„  (22') 

Jc=Cf+C, 

where  the  integral  is  understood  to  be  of  Cauchy’s  principal  value  type.  The  first  term  on 
the  right  can  be  regarded  as  the  physical  optics  approximation  to  the  surface  field  in  the 
illuminated  region. 

The  far-zone  (r  — >  00)  fields  can  be  determined  by  using  the  principal  asymptotic 
forms  for  the  various  Hankel  functions.  For  the  quantity  ^11,  far-zone  approximation  can 
be  obtained  by  deforming  the  path  of  the  integral  in  (33)  and  evaluating  by  the  saddle 
point  method.  The  result  is 


OT  (x,y;Xo,t/o) 


nr2  (Ao  + 


Having  provided  the  integral  representation  and  expressions  for  the  various  fields,  we 
present  in  the  next  section  details  on  the  MEI  method  of  terminating  the  computational 
domain  in  the  finite  scheme. 


Ild.  MEI  Method: 

As  already  remarked,  the  MEI  method  allows  one  to  generate  near-field  conditions  for 
the  scattered  field  V’  to  simulate  free  space  at  the  terminating  boundary.  We  will  describe 
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the  procedure  to  determine  the  coefficients  a*  in  (7).  Mei  et  al.,  [4]  postulated  that  there 
exists  a  linear  combination  of  the  scattered  field  over  a  small,  discrete,  spatial  domain,  pk, 
such  that 

K 

+  ak‘^k  =  0-  (45) 

The  coefficients  are  postulated  to  be  dependent  upon  the  location  of  the  field  and 
geometry  of  the  scatterer,  but,  independent  of  the  excitation  of  the  scattered  field.  The 
last  postulate  enables  one  to  determine  the  coefficients  ak  using  a  finite  set  of  linearly 
independent  scattered  fields  caused  by  different  excitations.  It  is  not  the  purpose  of  the 
present  report  to  test  the  validity  of  these  postulates.  Rather,  we  employ  this  method  to 
develop  a  relatively  fast  numerical  method  for  predicting  both  the  sky  wave  as  well  as  the 
ground  wave  for  propagation  over  inhomogeneous  and  irregular  terrain.  The  starting  point 
is  an  accurate  representation  of  near  fields  such  as  equation  (22).  The  scattered  field  is 
given  from  (22)  as 

v^(^  =  -  ffo.(/5)  =  -J  (46) 

c 


Making  use  of  this  in  (45)  we  see  that 


I  H,{p')  ^(p',po)  +  ^afc^(p',pit)  di'  =  0. 
•L  L  k-\ 


The  specific  excitations  (which  they  termed  as  metrons)  suggested  by  Mei  et  al,  were 


H,  =  l,  cos(27rs),  sin(27rs),  cos(47rs),  and  sin(47rs), 


(48) 


where  s  is  the  normalized  arc  length  that  varies  between  0  and  1  on  the  scatterer.  We  will 
label  these  as  sinusoidal  metrons.  We  were  not  able  to  consistently  generate  meaningful 
results  using  sinusoidal  metrons,  and  thus  considered  other  choices.  A  closer  look  at  (47) 
provides  insight  as  to  why  the  sinusoidal  metrons  cannot  be  expected  to  work  all  the  time. 
Since  the  coefficients  a*  have  been  postulated  to  be  independent  of  the  excitation  (which 
determines  on  the  scatterer),  it  is  implied  that  the  terms  within  the  bracket  be 

identically  zero  for  all  points  p'  on  the  scatterer.  The  term  g(p',pk)  is  the  field  at  p'  due 
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to  a  line  source  located  at  pk  ■  Hence  we  require  that  the  field  due  to  a  linear  combination 

of  the  these  line  sources  (K  +  1  in  number)  vanish  pointwise  on  the  scatterer.  But  this 

would  be  difficult  to  accomplish  as  there  are  only  a  few  discrete  sources  {K  is  typically  4-5 

depending  on  the  type  of  the  computational  molecule  chosen).  This  is  particularly  true 

for  a  large  scatterer.  A  more  reasonable  criterion  is  to  require  that  the  field  be  minimum 

in  some  sense  on  the  scatterer.  We  choose  the  coefficients  ajt  so  that  integrated  square 

residual  „ 

K  ^ 

R  =  I  g(P,M  +  J2akQ{P,Pk)  di'  (49) 

C 

is  a  minimum  on  the  scatterer.  In  other  words,  we  determine  the  coefficients  by  requiring 
that  =  0.  This  position  was  first  taken  by  Jetvic  and  Lee  [16],  who  considered  the 
case  of  perfectly  conducting  cylinder  to  demonstrate  the  success  of  this  approach.  This 
criterion  results  in 

I  S(?,hWG'J.)'U.'  =  -  f  Q{pMS\?Jn)dl\  n  =  1,2 . K.  (60) 

k=l  c  C 

The  coefficients  ajt  can  be  obtained  by  solving  the  linear  system  defined  by  (50).  We  shall 
label  the  particular  choice  of  metrons  in  (50)  as  the  Q*  metrons.  The  coefficient  matrix 
in  (50)  is  Hermitian  symmetric.  The  number  of  integrals  required  per  a  five-point  star 
mesh  {K  =  4)  is  14.  If  every  node  is  utilized  in  the  evaluation  of  an  integral,  the  total 
number  of  flops  required  to  perform  the  integration  per  a  five  star  mesh  is  14iV.  The 
number  of  Green’s  function  evaluations  per  a  five-point  star  mesh  is  bN .  The  total  time 
required  to  fill  the  coefficient  matrix  in  (50)  for  N  nodes  on  the  object  boundary  is  then 
lAN'^tf  -j-  bN^tg  =  where  tj  and  tg  are  the  respective  times  required  per  flop  and 

a  single  evaluation  of  the  Green’s  function.  For  example,  if  it  takes  1  millisecond  for  a 
single  evaluation  of  the  Green’s  function,  the  total  time  spent  in  the  evaluation  of  the 
Green’s  functions  for  N  =  1000  would  be  83  minutes.  It  is  therefore  very  important  to 
minimize  the  time  required  to  calculate  the  Green’s  function.  The  total  time  needed  to 
fill  the  matrices  (11)  through  (14)  is  still  of  order  0{N^).  The  fill  time  may  be  reduced 
approximately  by  a  factor  of  half  by  skipping  every  other  point  in  the  evaluation  of  the 
integrals.  The  next  section  deals  with  the  inversion  of  the  block  system  defined  in  (8-10). 
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lie.  Inversion  hy  Ricatti  Transform 

The  overall  matrix  generated  by  the  block  system  of  equations  in  (8-10)  is  highly 
sparse  and  is  almost  tridiagonal.  It  may  be  inverted  efficiently  by  the  Ricatti  transform 
method  [8],  which  was  originally  described  in  [17].  To  this  end  we  write 

+  S„,  n  =  2,...,N,  (51) 

where  the  transform  matrix  R„  is  of  order  M  x  M,  and  the  shift  vector  S„  is  of  order 
M  X  1.  Substituting  into  equation  (9),  we  get 

-  [A„R„  +  B„]-'  +  [A„R„  +  [F"  -  A„S„] .  (52) 

We  then  observe  by  comparing  with  (51)  that  for  n  =  2, . . . ,  N  —  1 

Rn+x  =  -  [A„R„  +  B„]-'  C„,  S„+i  -  [A„R„  +  B„]“'  [F"  -  A„S„] .  (53) 

The  end  equation  (8)  may  be  used  to  determine  the  transform  matrix  R2  and  S2. 
Substituting  (51)  into  (8)  and  carrying  out  some  algebraic  manipulations,  we  get 

R2  =  (AiC2-'A2  -  Bi)"'  (Cl  -  AiC2-'B2)  ,  (54) 

and 

S2  =  (A2  -  C2A-iBa)"'  (F^  -  C2A-'f1)  .  (55) 

The  higher  order  matrices  can  be  determined  from  the  above  two  by  using  the  recursive 
relations  in  (53).  Finally,  by  using  the  other  end  equation  (10)  in  conjunction  with  (51), 
we  solve  for  as 

=  [{An  +  CnRn-i)^n  +  Bat]"'  [F^  -  C^vSn-i  -  (A^  +  C^Riv-OSiv]  •  (56) 

The  remaining  vectors  can  then  be  determined  by  the  tranformation 

equations  in  (51).  Each  of  the  matrices  A„,  C„,  Sn,  and  has  M  non-zero  elements 
while  B„  has  3M  non-zero  elements.  The  matrix  R„  is  dense  and  has  elements.  The 
total  number  of  non-zero  elements,  that  need  to  be  stored  with  this  algorithm  is  then 

Ns  =  {N-  1){M^  +  M)  +  6MN  ~  (M  +  7)MN  =  0{N).  (57) 
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For  example,  with  N  =  1000  and  M  =  5,  the  algorithm  needs  around  0.96  MB  of  RAM  for 
implementation  in  double  precision.  If  N  is  now  increased  to  10,000,  the  required  memory 
is  increased  to  9.6  MB.  This  is  in  contrast  to  the  integral  equation  methods,  where,  a 
scaling  in  A"  by  a  factor  of  10  results  in  100  fold  increase  in  memory.  The  inversion  time 
of  the  algorithm  is  only  of  order  0(iV).  The  overall  CPU  time  of  the  present  method  is 
dominated  by  the  matrix  fill  time  discussed  in  the  previous  section. 
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III.  NUMERICAL  RESULTS 


The  entire  code  was  developed  in  double  precision  in  FORTRAN.  The  linear  equations 
generated  by  (4)  and  (5)  for  the  determination  of  the  coefficients  were  solved  by  Gaussian 
elimination.  Due  to  the  generation  of  a  highly  ill-conditioned  matrix,  Gaussian  elimination 
is,  however,  not  reliable  for  the  inversion  of  the  matrix  associated  with  the  coefficients  a*,. 
By  treating  it  as  a  linear  least  square  problem,  we  have  used  Singular  Value  Decompsition 
(SVD)  [19]  to  construct  a  minimum  norm  solution  to  the  matrix  equation.  SVD  is  also 
appropriate  for  the  case  of  sinusoidal  matrons  where  we  generate  more  equations  that  the 
number  of  unknowns.  The  effective  rank  of  the  matrix  is  determined  by  treating  as  zero 
those  singular  values  which  are  less  that  a  predetermined  number  ‘Rcond’  times  the  largest 
singular  value.  The  condition  number,  K2-,  of  the  matrix  in  the  2-norm  is  the  ratio  of  the 
highest  singular  value  to  the  lowest  singular  value.  Of  course,  if  Rcond  is  set  less  than  the 
infimum  of  1/k2  over  all  the  nodes  on  the  boundary,  SVD  uses  all  the  singular  vectors  and 
produces  the  same  result  as  Gaussian  elimination  for  a  square  matrix.  Integration  in  (46) 
was  performed  by  the  Simpson’s  rule  either  using  all  nodes  on  the  scatterer  or  using  every 
other  node.  The  latter  reduces  the  integration  time  by  a  factor  of  half.  The  geometry  of  the 
scatterer  was  specified  in  a  discrete  form  by  the  nodes  at  which  the  unknowns  are  defined. 
Interpolation  with  quadratic  elements  [20]  was  used  to  generate  a  continuous  object.  In 
this  way  the  code  was  capable  of  handling  a  rather  arbitrary  geometry.  The  various  Hankel 
functions  with  complex  arguments  in  (37)  and  (39)  were  generated  by  implementing  the 
algorithm  of  du  Toit  [15]. 

To  verify  calculations  by  the  complex  image  approach  of  section  lie,  we  first  present 
results  on  the  computation  of  the  incident  magnetic  field.  Fig.  6  shows  the  magnitude  of 
the  normalized  magnetic  field  due  to  a  vertically  polarized  line  source  on  the  surface  of  a 
flat,  lossy  plane.  The  ground  constants  correspond  to  Cr  =  15  and  C7r  =  6.  The  magnetic 
field  is  normalized  to  the  free-space  value,  which  is  true  of  all  the  results  shown  in  the 
report.  The  constant  To  was  chosen  to  be  at  least  5  to  make  the  approximation  in  (37) 
valid.  Fourteen  complex  images  were  chosen  to  approximate  the  spectral  function,  although 
10  were  found  to  produce  identical  results.  Expression  (37)  was  used  for  r’2  <  50,  while 
(39)  was  used  otherwise.  The  upper  limit  for  u  on  the  steepest  descent  path  in  the  integral 
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(39)  is  determined  such  that  the  imaginary  part  of  the  complex  distance  D  takes  a  value 
of  5.  The  arc  length  value  in  Fig.  6  is  equal  to  zero  right  under  the  source.  Comparision  is 
shown  with  the  calculation  based  on  Sommerfeld  integral  given  in  [1].  Excellent  agreement 
is  seen  between  the  two,  thus  validating  the  complex  image  aproach.  Several  other  cases 
have  also  been  successfully  tested  including  the  cases  of  large  conductivity  such  as  =  10, 
and  Gr  =  200.  Where  applicable,  we  used  Nt  =  10  and  To  =  5  in  the  numerical  results 
shown  in  this  report. 

We  next  compare  the  numerical  results  for  propagation  over  a  semicircular  boss  on 
a  perfectly  conducting  plane  (Aq  =  0),  for  which  the  exact  solution  in  terms  on  cylin¬ 
drical  harmonics  exists  [1].  Results  are  shown  both  for  the  sinusoidal  metrons  and  the 
G*  metrons.  A  vertically  polarized  source  (9o  =  90)  was  assumed  to  be  located  at 
(-0.75A,0.1A)  near  a  semicircular  boss  of  radius  0.5A  with  center  at  the  origin.  Six  layers 
(M  =  4,  recalling  that  M  denotes  the  number  of  interior  layers  not  counting  the  object 
and  outer  boundaries)  with  an  inter-layer  spacing,  h,  of  0.05A  were  used  to  discretize  the 
computational  domain.  The  distance  between  the  outermost  layer  and  the  object  bound¬ 
ary  is  (M  +  l)h  =  0.25A.  The  outer  boundary  was  discretized  into  roughly  twenty  segments 
per  wavelength  resulting  in  =  49.  Rcond  for  SVD  was  set  at  l.D-6.  The  arc  length,  S, 
takes  a  value  0  at  the  left  end  of  the  boss  and  increases  in  the  clockwise  direction.  Fig. 
7  shows  a  good  agreement  for  the  magnitude  of  the  surface  magnetic  fields  between  the 
numerical  and  the  exact  results  both  for  the  sinusoidal  and  G*  metrons.  The  supremum, 
(indicated  as  ‘Con’  in  the  plots),  of  the  2-norm  condition  number  of  the  4  X  4  coefficient 
matrix  in  the  case  of  G*  metrons,  and  the  5  x  4  matrix  in  the  case  of  sinusiodal  metrons 
over  all  nodes  on  the  boundary  is  also  indicated  on  the  plot.  Fig.  8  shows  a  comparison 
of  the  relative  pattern  of  the  source  in  the  presence  of  the  boss.  Good  agreement  is  seen 
except  near  the  grazing  angles  in  the  shadow  region  where  a  discrepancy  of  about  2  dB 
is  seen.  The  accuracy  of  the  solution  can  be  affected  by  varying  Rcond,  M,N  and  h. 
While  it  is  generally  true  that  the  accuracy  improves  when  M  or  N  is  increased  or  when 
is  h  decreased,  difficulty  to  compute  the  various  quantities  precisely  over  a  small,  discrete, 
spatial  domain  tends  to  obscure  this.  The  last  two  parameters  affect  the  discretization 
errors,  while  M  and  h  together  influence  the  distance  between  the  object  boundary  and 
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the  outer  boundary.  All  other  paramaters  being  fixed,  if  M  is  increased  and  h  is  decreased 
while  maintaining  a  constant  seperation  between  the  object  and  outer  boundaries,  the  so¬ 
lution  improves  only  slightly  as  indicated  in  Figs.  9  and  10.  In  Fig.  9  we  use  M  =  6  and 
h  =  0.0357A,  while  in  Fig.  10  we  use  M  =  8  and  h  =  0.0278A.  The  slight  improvement 
over  the  results  of  Fig.  7  is  thought  to  be  due  to  a  decreased  discretization  error  achieved 
by  the  use  of  smaller  h.  Compared  to  Fig.  7  it  is  also  seen  that  increases  in  both  cases. 
This  is  due  to  the  larger  size  of  the  matrices  involved.  The  effect  of  integration  on  the 
numerical  solution  was  also  investigated.  In  one  case,  integration  was  performed  by  using 
all  nodes  on  the  boundary,  while  in  the  other  case  every  other  node  was  used.  The  former 
is  labeled  as  full  integration,  while  the  latter  as  half  integration.  Fig.  11  shows  the  results 
with  full  and  half  integration.  The  solution  is  less  sensitive  to  integration  for  the  G*  case 
than  it  is  for  the  sinusoidal  metron  case.  This  trend  has  been  noticed  for  other  shapes 
as  well.  To  see  the  effect  of  node  density  on  the  solution,  the  number  of  segments  on  the 
outer  boundary  was  increased  to  29  per  wavelength.  This  was  with  a  view  to  increase 
the  accuracy  of  the  solution.  Fig.  12  shows  the  solution  with  the  higher  node  density. 
The  G*  metron  solution  appears  to  have  improved  slightly,  whereas  the  sinusoidal  metron 
case  deteriorated  slightly  when  compared  with  the  20  segment /wavelength  case  of  Fig.  9. 
Note  the  increase  Q,  for  both  cases  when  compared  to  Fig.  9.  This  is  because  the  fields 
calculated  over  a  small  spatial  domain  tend  to  be  more  linearly  dependent  than  over  a 
larger  domain. 

Parametric  studies  were  also  made  for  larger  cylinders.  Fig.  13  shows  the  results 
for  a  radius  of  5A.  With  the  indicated  choice  of  parameters,  the  number  of  nodes  on  the 
boundary  is  N  =  331.  Once  again  the  agreement  between  numerical  and  exact  results  is 
good.  Some  spurious  oscillations  are,  however,  seen  in  the  numerical  solution  in  the  last 
part  of  the  shadow  region.  The  total  computational  time  needed  on  a  80486-50  Personal 
Computer  for  the  calculation  of  ground  wave  data  at  331  points  and  the  far-zone  data  at 
181  angles  utilizing  every  other  point  for  integartion  was  2|  minutes.  Almost  all  of  the 
time  was  required  for  filling  the  various  matrices,  with  the  inversion  taking  only  a  paltry 
1  second.  Fig.  14  shows  the  relative  pattern  of  the  source.  It  is  seen  that  a  slightly  better 
agreement  is  obtained  with  the  sinusoidal  metrons  than  with  the  G*  metrons.  As  in  the 


21 


previous  case  the  accuracy  improved  only  slightly  when  the  number  of  interior  layers  is 
increased  without  changing  the  distance  between  the  object  and  outer  boundaries.  This  is 
clearly  seen  from  Figs.  15  and  16  which  show  results  for  M  =  6  and  M  —  8  respectively. 
For  the  sinusoidal  case,  however,  the  result  with  M  =  8  is  once  again  less  accurate  when 
compared  with  M  =  4.  Increasing  Rcond  had  no  effect  on  the  solution.  However,  when 
the  number  of  segments  on  the  outer  boundary  is  increased  to  34  per  wavelength,  a  more 
accurate  solution  is  obtained  as  seen  in  Fig.  17.  At  this  node  density  the  spacing  between 
nodes  on  the  object  boundary  is  the  same  as  h.  Fig.  18  shows  the  numerical  solution 
with  the  higher  node  density  of  34  segments  per  wavelength  for  M  =  8  and  h  =  0.0278A. 
For  the  G*  metron  case  the  solution  remains  practically  unchanged  from  that  of  Fig. 
16,  but  the  result  for  sinusoidal  metron  case  is  more  meaningful  with  the  higher  node 
density.  The  effect  of  moving  the  terminating  boundary  further  away  from  the  object 
boundary  is  seen  by  examining  Fig.  19.  For  the  parameters  indicated,  the  outer  boundary 
is  at  a  distance  of  0.45A  from  the  cylinder.  For  the  G*  metron  case,  the  hump  in  the 
solution  around  S'  =  70  is  still  present,  although  an  improvement  is  seen  at  the  left  end 
of  the  boss  (the  illuminated  side).  Overall,  the  G*  solution  agrees  better  with  the  exact 
solution  when  compared  with  Fig.  13  where  the  seperation  between  the  boudaries  is  0.25A. 
Increasing  the  node  density  to  34  per  wavelength  did  not  improve  the  solution  any  further. 
Further  improvement  in  the  solution  can  be  affected  by  decreasing  h  which  decreases  the 
discretization  errors.  This  is  clearly  seen  from  Fig.  20  where  the  solution  with  h  =  0.035A 
compares  better  than  the  one  with  h  =  0.05A.  The  spurious  oscillations  can  be  reduced 
by  (i)  increasing  the  seperation  between  the  object  and  the  terminating  boundaries,  and 
(ii)  concurrently  discarding  those  singular  vectors  that  are  corrupted  by  roundoff  errors. 
The  latter  is  accomplished  by  decreasing  Rcond  during  the  SVD  solution  of  the  coefficients 
ak-  In  general,  it  has  been  found  that  for  a  given  spatial  resolution,  the  coefficient  matrix 
becomes  more  ill-conditioned  (as  evidenced  by  a  higher  condition  number)  as  the  seperation 
between  the  object  and  outer  boundaries  gets  larger.  Consequently,  Rcond  has  to  be  set 
a  lower  value  in  order  to  achieve  a  meaningful  solution.  Fig.  21  shows  the  solution  for  a 
separation  of  0.6A  and  Rcond=l.D-4.  The  accuracy  is  seen  to  be  the  best  of  all  the  results 
shown  so  far.  If  however,  only  an  approximate  solution  is  desired,  the  parameters  of  Fig. 
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13  are  sufficient.  The  results  seem  to  suggest  that  the  accuracy  depends  primarily  on 
the  distance  between  the  object  boundary  and  the  outer  bormdary  provided  that  various 
quantities  are  calculated  precisely.  Twenty  segments  per  wavelength,  six  layers,  and  an 
interlayer  spacing  of  A/20  seem  to  produce  reasonable  results  for  most  cases.  If  G*  metrons 
are  used,  integration  can  be  done  using  every  other  point. 

Next,  we  consider  an  obstacle  that  has  both  concave  and  convex  potions.  A  good 
candidate  is  the  case  of  a  gaussian  hill  for  which  results  are  available  in  the  literature  [18]. 
Such  a  hill  has  also  been  considered  in  [1].  The  height  of  the  hill,  hg{x),  over  the  y  =  0 
plane  is  given  by 

hgix)  = 

We  have  choosen  A  =  C  -  lOA/3,  5  =  0.  The  obstacle  is  taken  as  the  portion  of  the 
gaussian  hill  where  the  height  exceeds  A/100.  For  the  chosen  parameters,  this  occurs  for 
lx|  ~  8A.  The  total  arc  length  of  the  obstacle  is  about  18A.  The  constitutive  parameters 
of  the  earth  are  taken  as  =  10,  <7  =  lOmS/m  at  a  frequency  of  IMHz.  The  latter 
corresponds  to  Or  =  180.  A  vertically  polarized  source  is  located  at  (—50A/3,  A/100)  to 
the  left  of  the  hill.  Fig.  22  shows  the  incident  magnetic  field,  Hoz,  on  the  hill  normalized 
to  the  free-space  value.  The  dashed  line  corresponds  to  the  field  over  a  flat  surface.  The 
quantity  ^11  in  (42)  has  been  calculated  using  (39)  since  r2  >  50  on  the  obstacle.  As 
expected,  the  field  matches  at  the  two  ends  of  the  obstacle  with  the  field  over  a  flat 
surface.  In  [18]  the  solution  is  obtained  via  Volterra  type  integral  equation  starting  from 
parabolic  approximation  of  the  Helmholtz  equation.  By  its  very  nature,  the  method  of  [18] 
considers  one-way  propagation  and  ignores  backscattering.  The  present  method,  on  the 
other  hand,  makes  no  such  approximations.  Fig.  23  compares  the  numerical  solution  of  the 
present  method  with  that  of  [18].  Both  the  G*  and  the  sinusoidal  metron  solutions  have 
been  obtained  with  20  segments  per  wavelength  node  density,  six  layers,  h  =  0.05A,  and 
Rcond=l.D-6.  For  these  values  of  parameters,  the  munber  of  imknowns  on  the  obstacle  is 
361.  Full  integration  was  employed  in  both  cases.  The  solution  obtained  with  G*  metrons 
agrees  fairly  well  the  results  of  [18],  whereas  the  one  with  sinusoidal  metrons  is  highly 
erroneous.  The  increased  field  strength  on  the  illuminated  side  is  due  to  focussing  by 
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the  concave  portion  of  the  hill.  The  G*  metron  solution  correctly  predicts  it.  No  such 
focussing  is  present  in  the  incident  field  as  can  be  observed  from  Fig.  22.  Thus,  physical 
optics  cannot  be  expected  to  predict  the  correct  field  on  the  surface.  It  is  to  be  noted 
that  the  erroneous  result  with  sinusoidal  metrons  is  not  due  to  round  off  errors,  for  the 
G*  solution  has  a  higher  condition  number  and  still  produces  a  good  result.  This  was  also 
checked  by  repeating  the  numerical  solution  with  various  values  of  Rcond.  The  sinusoidal 
metron  solution  remained  erroneous  irrespective  of  the  value  of  Rcond.  In  an  attempt  to 
improve  the  solution,  a  higher  node  density  of  30  segments  per  wavelength  was  also  tried. 
However,  the  solution  only  worsened.  It  is  concluded  that  the  sinusoidal  metrons  are  not 
satisfactory  for  arbitrary  shapes  and  that  they  do  not  lead  to  the  correct  values  of  the  MEI 
coefficients.  Choice  of  metrons  is  not  arbitrary  as  the  originators  of  the  MEI  method  claim. 
The  G*  metrons  have  a  physical  significance  in  that  they  represent  the  fields  generated  by 
line  sources  in  an  environment  compatible  with  the  original  problem.  No  such  statement 
can  be  mode  concerning  the  sinusoidal  metrons. 

It  is  interesting  to  compare  the  merits  of  the  present  method  with  G*  metrons  vis  — 
a  -  vis  the  method  of  [1]  which  solves  the  problem  using  a  magnetic  field  integral  equation 
and  boundary  element  method.  It  took  5  minutes  and  20  seconds  on  a  80486-50  PC  to 
produce  the  results  for  the  ground  wave  data  shown  in  Fig.  5  as  well  as  to  generate  sky 
wave  data  at  181  angles.  Identical  results  were  obtained  with  half  integration  which  took 
only  3  minutes.  This  is  in  contrast  to  the  method  of  [1]  which  took  75  minutes  for  the  same 
problem.  Of  course  [1]  uses  the  free-space  Green’s  function  in  contrast  to  the  half-space 
Green’s  function  employed  here.  The  number  of  unknowns  in  [1]  was  725  (obtained  with 
12  nodes  per  wavelength).  For  the  same  number  of  unknowns,  the  present  method,  being 
of  order  would  take  about  21  minutes,  which  is  still  faster  by  a  factor  of  about 

four.  If  half  integration  is  used  the  method  is  faster  by  a  factor  of  eight.  The  savings  in 
memory  in  the  present  method  are  tremendous.  For  the  problem  at  hand,  the  method  of 
[1]  with  725  unknowns  requires  a  storage  space  of  at  least  8.4  MB  compared  to  only  254 
kB  needed  with  the  present  method.  To  speed  up  the  calculation  of  the  Green  s  function 
we  have  assumed  =  0  in  (43).  This  amounts  to  assuming  that  the  terrain  outside  the 
obstacle  is  made  up  of  a  perfect  conductor.  This  results  in  substantial  savings  in  time. 
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For  example,  it  took  3  hours  and  10  minutes  on  a  80486-50  PC  to  solve  the  problem  of 
gaussian  hill  using  the  exact  Green’s  function  and  half  integration,  compared  to  only  3 
minutes  with  the  approximate  Green’s  function!  The  results  were  identical  except  near 
the  right  end  of  the  hill  where  the  two  differed  slightly.  The  value  at  x  —  Xq  =  160  was 
1.2  with  the  exact  Green’s  function  compared  to  about  0.9  with  the  approximate  Green’s 
function.  One  would  expect  the  use  of  approximate  Green’s  function  to  be  the  worst  when 
Ao  is  significantly  different  from  zero.  We  considered  the  case  of  highly  lossy  earth  with 
Cj.  =  5  and  a  =  0.056mS/m  at  1  MHz  =  1).  A  gaussian  hill  with  A  =  X,  B  =  0,  and 
C  =  0.76A  was  considered.  A  vertically  polarized  source  is  placed  at  (— 3A,A/10).  Six 
layers  were  chosen  with  h  =  0.05A  and  integartion  was  performed  utilyzing  all  the  nodes. 
Fig.  24  shows  the  calculations  with  the  exact  and  the  approximate  Green’s  functions.  The 
field  that  would  exist  on  flat  earth  is  also  shown.  Calculations  with  the  exact  Green’s 
function  took  30  seconds  while  those  with  the  exact  Green’s  function  took  45  minutes. 
However,  it  is  seen  that  the  results  do  not  differ  much  from  each  other,  justifying  the  use 
of  the  approximate  Green’s  function. 

It  is  possible  to  increase  the  accuracy  of  the  numerical  results  shown  in  Fig.  23  by 
increasing  the  seperation  between  the  obstacle  and  the  outer  boundary.  Fig.  25  shows 
calculations  with  6,  8,  and  10  layers  all  with  h  =  0.05A.  It  is  seen  that  the  result  with 
10  layers  is  generally  in  best  agreement  with  the  results  of  [18].  The  numbers  indicated 
within  the  paranthesis  in  the  caption  are  the  values  of  Rcond  assigned  in  each  case  to  the 
first  15  wavelengths  and  the  last  3  wavelengths  respectively  on  the  boundary.  As  discussed 
previously,  the  condition  number  increases  generally  as  the  seperation  is  increased.  We 
selectively  filter  out  the  singular  vectors  which  are  highly  corrupted  by  roundoff  errors  by 
choosing  a  higher  value  of  Rcond.  The  far-zone  fields  are  however  insensitive  to  the  slight 
variations  of  the  surface  fields.  Fig.  26  shows  the  relative  pattern  of  the  source  with  6  and 
8  layers.  It  is  seen  that  the  two  cases  produce  virtually  the  same  results.  The  radiation 
pattern  of  the  source  in  the  direction  of  the  hill  is  highly  perturbed  by  the  presence  of 
the  latter.  Away  from  the  hill  the  relative  pattern  in  the  presence  of  the  hill  is  not  very 
different  from  the  pattern  of  a  vertical  source  in  free-space  except  near  grazing  angles 
where  it  has  a  deep  minimum. 
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Finally,  we  present  results  on  the  effect  of  relative  location  of  the  source  with  respect 
to  the  obstacle  on  the  fields.  A  vertically  polarized  source  was  placed  near  a  gaussian  hill 
having  A  =  =  0,  and  C  =  A.  The  total  arc  length  of  the  obstacle  is  5.1A.  The 

ground  constants  were  e,  =  14.9,  a.  =  14.58.  Six  layers  with  20  segments  per  wavelength, 
and  h  =  0.05A  were  chosen  in  the  numerical  model.  In  one  case  the  source  was  placed  at 
the  bottom  of  the  hill  at  a  horizontal  distance  of  4.53A  from  the  peak.  In  the  second  case 
the  source  was  placed  at  the  top  of  the  hill  (xo  =  0).  Fig.  27  shows  the  normalized  ground 
wave  on  the  hill.  Once  again  there  is  focussing  on  the  illuminated  side  of  the  hill  when 
the  source  is  at  the  bottom.  The  field  for  the  source  at  the  top  of  the  hill  is  symmetric 
as  expected.  Note  that  the  field  in  the  shadow  region  of  the  hill  for  lo  =  -4.53A  is  as 
strong  as  the  corresponding  field  for  Xo  =  0.  There  is  no  apparent  advantage  of  siting  the 
antenna  at  the  top  of  the  hill  to  receive  the  ground  wave.  However,  the  radiation  pattern 
of  the  source  from  (f)  =  0  to  90°  is  significantly  affected  by  the  presence  of  the  hill.  This 
can  be  seen  from  Fig.  28  which  shows  the  relative  pattern  for  the  two  source  locations. 
Irrespective  of  the  source  location,  the  radiation  pattern  has  a  deep  minimum  near  grazing 
angles  as  is  characteristic  of  vertical  antennas  over  lossy  earth. 
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IV.  SUMMARY 


A  fast,  finite  difference  method  that  includes  all  aspects  of  wave  phenomenon  such 
as  reflection,  refraction,  diffraction,  and  backscattering  is  presented  for  predicting  two  di¬ 
mensional  propagation  over  inhomogeneous,  irregular  terrain.  The  terrain  is  characterized 
by  its  elevation  and  impedance,  which,  in  turn  depends  on  the  ground  constants  of  the 
earth.  Both  of  these  may  vary  with  distance.  The  terrain  topology  data  is  specified  at 
discrete  points.  Interpolation  using  quadratic  elements  is  done  to  define  a  continuous  ge¬ 
ometry.  The  computational  domain  for  the  problem  consists  of  the  area  bounded  by  the 
terrain  at  the  bottom  and  a  truncating  boundary  at  the  top.  To  simulate  free-space  on  the 
truncating  boundary,  discrete,  near-field  radiation  condition  of  Mei  type,  derivable  from 
an  integral  representation  of  the  fields,  is  imposed.  Green’s  function  for  half-space  is  used 
to  reduce  the  number  of  unknowns.  Unknowns  are  distributed  on  the  terrain  only  where 
its  elevation  is  non-zero  and/or  where  its  impedance  differs  from  a  reference  value.  The 
computational  domain  is  discretized  using  interior  layers  between  the  truncating  and  the 
object  boundaries.  Finite  difference  coefficients  valid  for  an  irregular,  non-othogonal  mesh 
are  presented.  Accurate  expressions  are  provided  for  the  Green’s  function  and  incident 
fields  over  a  constant-impedance,  flat  plane.  The  expressions  permit  rapid  computation  as 
they  do  not  involve  the  troublesome  Sommerfield  integrals.  Results  are  presented  for  the 
ground  wave  as  well  as  the  sky  wave. 

The  truncating  boundary  can  be  in  the  near-field  of  the  obstacle;  as  near  as  a  A/4 
away  from  it,  and  the  method  works  both  for  concave  and  convex  geometries.  Good  results 
have  been  obtained  with  a  node  density  of  about  20  per  wavelength  and  an  interlayer 
spacing  of  about  A/20.  The  metrons  used  for  the  determination  of  the  Mei  coefficients  are 
proportional  to  the  complex  conjugate  of  the  Green’s  function  and  fully  accomodate  the 
environment  of  the  problem.  Singular  value  decomposition  which  permits  filtering  of  the 
space  spanned  by  the  corrupted  singular  vectors  is  used  to  solve  for  the  coefficients. 

The  method  is  attractive  for  large  terrain  obstacles  where  other  methods  tend  to 
be  slow.  The  overall  computational  time  of  the  method  is  of  order  O(iV^).  Storage  re¬ 
quirements  are  of  order  0(V),  where  N  is  the  total  number  of  unknowns  (nodes).  As  an 
example,  for  a  terrain  obstacle  extending  over  18  wavelengths,  the  number  of  nodes  on  the 
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boundary  is  361.  To  compute  ground  wave  data  at  all  of  these  nodes  as  well  as  to  compute 
the  far-fields  at  181  angles,  the  method  takes  around  3  nainutes  on  a  80486-50  PC.  The 
total  storage  required  for  the  problem  is  around  254  kB.  It  would  take  around  300  minutes 
with  a  storage  space  of  around  2.5  MB  to  solve  the  problem  of  a  terrain  where  the  obstacle 
extends  over  180  wavelengths.  At  a  frequency  of  1  MHz  the  terrain  obstacle  corresponds 
to  one  whose  arc  length  is  54  km.  At  a  frequency  of  1  GHz  it  corresponds  to  an  arc  length 
of  54  m. 

Certain  spurious  oscillations  and  over  estimation  of  the  field  have  been  observed  in 
some  regions  of  the  shadow  region.  This  is  attributed  to  the  fact  that  it  is  difficult  to 
calculate  near-fields  precisely  over  an  electrically  small,  discrete  domain.  These  extraneous 
effects  can  be  some  what  reduced  by  increasing  the  seperation  between  the  obstacle  and 
the  truncating  boundary.  It  is  speculated  that  they  can  be  further  reduced  by  choosing  a 
computational  molecule  that  accomodates  higher  order  derivatives.  More  efficient  schemes 
of  evaluating  the  Green’s  function  and  the  near-fields  will  have  to  be  made  to  further 
reduce  the  computational  time. 

Extension  to  the  three  dimensional  case  is  trivial  in  principle,  although  not  so  com¬ 
putationally.  Savings  over  integral  equation  methods  will  be  even  more  dramatic  in  the 
case  of  three  dimensional  obstacles.  Although  the  Green’s  function  using  complex  images 
will  be  little  more  involved  than  for  the  two-dimensional  case,  it  will  have  the  computa¬ 
tional  advantage  of  being  expressible  in  terms  of  an  exponential  functions  instead  of  Hankel 
functions  with  complex  arguments.  The  fields  decay  more  rapidly  with  distance  in  three 
dimensions  than  in  two  dimensions.  As  a  result  integration  could  be  done  little  more  effi¬ 
ciently  in  the  former  case.  The  corresponding  computational  molecule  will  now  consist  of 
seven  nodes  instead  of  five  in  two  dimensions.  Finite  difference  and  terminating  boundary 
conditions  will  have  to  be  developed  for  vector  fields  instead  of  scalar  fields.  These  task 
are  presently  being  pursued. 
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Surface  Magnetic  Field  on  Circular  Boss 
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Source  Pattern  Near  Circular  Boss 
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Surface  Magnetic  Field  on  Circular  Boss 
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Surface  Magnetic  Field  on  Circular  Boss 
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Surface  Magnetic  Field  on  Circular  Boss 
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Surface  Magnetic  Field  on  Circular  Boss 
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Source  Pattern  Near  Circular  Boss 


FIG.  14. 


44 


Surface  Magnetic  Field  on  Circular  Boss 
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Surface  Magnetic  Field  on  Circular  Boss 
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Surface  Magnetic  Field  on  Circular  Boss 
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Surface  Magnetic  Field  on  Circular  Boss 
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Surface  Magnetic  Field  on  Circular  Boss 
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Surface  Magnetic  Field  on  Circular  Boss 
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Surface  Magnetic  Field  on  Gaussian  Hill 
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Surface  Magnetic  Field  on  Gaussian  Hill 
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Surface  Magnetic  Field  on  Gaussian  Hill 
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Pattern  of  Vertical  Source  Near  Gaussian  Hill 
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Pattern  of  Vertical  Source  Near  Gaussian  Hill 
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